/**
 * 计算数字滤波器的频率响应
 * num 是数字滤波器的分子多项式系数
 * den 是数字滤波器的分母多项式系数
 * num_order 是分子多项式的阶数
 * den_order 是分母多项式的阶数
 * sign = 0 时，x_out 为频率响应的实部， y_out 为频率响应的虚部
 * sign = 1 时，x_out 为频率响应的模， y_out 为频率响应的幅角
 * sign = 2 时，x_out 为以 dB 为单位的频率响应， y_out 为频率响应的幅角
 * len 为，频率响应的取样点数
 */
void gain(double num[], double den[], int num_order, int den_order,
          double x_out[], double y_out[], int len, int sign)
{
    int i, k;
    double zr, zi;
    double re, im;
    double ar, ai, br, bi;
    double numr, numi;
    double de, temp;
    double freq;
    for(k = 0; k < len; k++)
    {
        freq = 0.5 * k / (len - 1);
        zr = cos(-2.0 * M_PI * freq);
        zi = sin(-2.0 * M_PI * freq);
        br = 0.0;
        bi = 0.0;
        for(i = num_order; i > 0; i--)
        {
            re = br;
            im = bi;
            br = (re + num[i]) * zr - im * zi;
            bi = (re + num[i]) * zi + im * zr;
        }
        ar = 0.0;
        ai = 0.0;
        for(i = den_order; i > 0; i--)
        {
            re = ar;
            im = ai;
            ar = (re + den[i]) * zr - im * zi;
            ai = (re + den[i]) * zi + im * zr;
        }
        br = br + num[0];
        ar = ar + 1.0;
        numr = ar * br + ai * bi;
        numi = ar * bi - ai * br;
        de = ar * ar + ai * ai;
        x_out[k] = numr / de;
        y_out[k] = numi / de;
        switch(sign)
        {
        case 1:
            temp = hypot(x_out[k], y_out[k]);
            y_out[k] = atan2(y_out[k], x_out[k]);
            x_out[k] = temp;
            break;
        case 2:
            temp = x_out[k] * x_out[k] + y_out[k] * y_out[k];
            y_out[k] = atan2(y_out[k], x_out[k]);
            x_out[k] = 10 *log10(temp);
            break;
        default:
            break;
        }
    }
}